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Mathematical modelling has emerged as a key process analysis tool in the area 
of Metals and Materials processing operations. Considerable advances have been made through 
mathematical modelling in terms of improved process efficiency and better product quality. 
Most of the models applied have been mechanistic in nature and are based on the conservation 
of heat, mass and momentum. Transport of heat mass and momentum in metal processing 
operations are characteristically multidimensional ( most often three dimensional ) and multi- 
phase ( slag, metal and gas ) in nature. Consequently, for effective simulation, three 
dimensional transport models (steady or transient) are required . The present work primarily 
concerns vvith the development of a steady state, two phase, three dimensional, turbulent flow 
calculation procedure from an existing two phase, two dimensional turbulent flow model. 

Thus, as a first step, the adequacy of an existing two dimensional turbulent flow model 
has been examined. This has been done by investigating melting phenomena in axisymmetrical 
gas stirred vessels for both aqueous and high temperature systems. Two dimensional 
predictions have been assessed against experimental measurements reported in literature and 
excellent agreement between the two demonstrated. Subsequent to this, the two dimensional 
calculation procedure has been upgraded to a steady, three dimensional configuration using the 
control volume based finite difference technique. The results from the three dimensional model 
were evaluated against standard benchmark solutions ( e.g., flow in a cubic cavity, entrance 
length against Reynolds number in ducts etc.) and physically realistic solutions were obtained 
for a wide range of conditions that are consistent with the observations reported in the literature. 
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Finally, the three dimensional turbulent flow model was used to predict the flow 
behaviour and RTD (residence time distribution) in a simple rectangular shaped tundish with 
no flow control devices. Predicted flows and RTD were qualitatively assessed against 
predictions reported in the literature and it was concluded that present model is capable of 
realistically simulating the hydro-dynamics of metallurgical tundish systems. 
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CHAPTER -1 


INTRODUCTION 


1.1 Introduction to the Thesis 

Modelling has played a key role in developing a holistic approach towards the 
understanding of metals and materials processing operations. Most of the metallurgical 
processes are characterised by high temperature, opacity of melts and large reactor size. These 
as a consequence, make the industrial scale metal processing units less than convenient case 
studies. Consequently, it has been customary to investigate the process dynamics of steelmaking 
with the aid of physical and mathematical models. 

Traditionally, physical modelling was the preffered tool [1] in metallurgy but of late, 
mathematical modelling has become more popular because of availability of several software 
packages, cheaper computational hardware, growing experience with the tackling of a broad 
range of computational problems and better understanding of the underlying fundamental 
concepts. With the development of new materials i.e., ceramics, composites etc. and a more 
challenging marketplace, the focus of a process metallurgist has undergone substantial 
reorientation and has shifted to cost control and better product quality. These warrant better 
process control and process optimisation . Towards this, mathematical modelling, over the last 
few years has played a key role, as both process control and process optimisation require the 
representation and understanding of the process / operation in quantitative form. 

Apart from lower costs, mathematical modelling has other inherent advantages over 
conventional forms of investigation viz. pilot plant scale experimentation and physical 
modelling. They are as follows: 



i. Speed : With the help of a mathematical model a investigator can study various 
configurations with ease, at a remarkable speed . 

ii. Information: The data derived from mathematical model is detailed and complete. The 
values of all possible variables in the domain of interest can be calculated and this 
information can be put to much more effective use. Mathematical modelling scores over 
other forms of investigation viz. physical modelling and pilot scale experimentation on 
this aspect. 

iii. Simulation flexibility: With the help of a mathematical model, small scale or cold 
models can be studied. Further, various theoretical possibilities such as, constant 
density and infinite reaction rate options can also be evaluated which at times help the 
researcher in developing a better understanding of the subject. 

The general methodology of a mathematical model development is shown 
schematically [2] in F/^.(l.l). This involves identifying the problem, carrying out simple 
calculations and then drawing a plan for detailed experimental and computational work. It 
needs to be stressed here that mathematical modelling together with physical modelling and 
pilot plant investigation needs to be carried out in a complementary fashion. As shown in 
F/g.(l.l) several iterations between modelling and experimental work may be required to 
develop a reliable mathematical model. 

Mathematical modelling, from a humble beginning has today developed into a very 
powerful tool capable of influencing our daily lives. Some of the notable successes of 
mathematical modelling have been in the accurate prediction of the election results based on 
extremely small samples, in weather forecasting and in the area of cold nuclear explosions. In 
the area of extractive metallurgy too, mathematical modelling has produced impressive results 
particularly in the last two decades. Some of the successful modelling examples are illustrated 
below : 

i. Aluminium extraction involves electrolytic deposition of alumina, dissolved in a 
molten salt using a graphite anode and a graphite cathode. A Hall-cell operation 
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F/g.(l.l) : The general methodology of mathematical model development [2]. 
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involves complex transport processes. The flow is unsteady and three dimensional and 
involves mass transfer in a multiphase ( metal, bath and gas ) system. Optimisation 
models [3,4] developed have not only helped in the evolvement of better cell design 
but also in the better understanding of the process which has resulted in vast 
improvements in the area of cell productivity and current efficiency [5]. 

ii. Continuous casting of steel is fast replacing the traditional ingot casting route as it is 
economically more attractive and also simplifies the process route by doing away with 
the maintainance intensive, primary rolling mills. However, the process of continuous 
casting involves solidification of the liquified metal in an oscillating mould which 
requires optimum process control. In the initial years, lack of proper understanding 
resulted in very high breakdowns and poor quality . Mathematical modelling of 
continuous casting helped in developing a better understanding of solidification in the 
mould and in designing of superior cooling systems for heat extraction from the mould 
[6,7]. This has resulted in controlling breakouts and in improving the quality of the 
continuous cast product. 

iii. Flash smelting of sulphide copper ores has inherent advantages of providing very high 
surface to volume ratio of the fine ore particles resulting in very high reaction rates. 
This has helped it in replacing the classical unwieldy reverberatory furnaces which 
provided very low surface area per unit volume and consequently, significantly lower 
processing rates. However, flash smelting is a very complex process involving a 
reacting two-phase, gas -solid mixture with several components. Mathematical 
modelling of this complex system helped in getting accurate information regarding 
temperature and concentration profiles [8,9] and of local reaction rates which helped in 
the design of larger and efficient flash smelters. 


1.2 Scope of the Present Work 

The brief introduction to the thesis indicates that mathematical modelling in recent years 
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has emerged as a major process analysis tool for various metal processing operations. Despite 
numerous modelling exercises taken up over the years [2-9], there remains considerable scope 
in broadening the applicability of the mathematical ( better called 'Process') models and their 
performance. A survey of the relevant literature indicates that mathematical modelling in the 
arena of metals processing has been most exhaustive in 2-D or axi-symmetrical systems. 
Transient two phase, 3-D calculations have been few and far between. Since, most metal 
processing operations do not involve symmetry, consequently the need for a three dimensional 
transport model can be readily visualized ( analysis of asymmetric gas injection, dual plug 
bubbling, flow in a tundish etc. necessitates a three dimensional turbulent flow calculation 
procedure). Similarly, although flow computations in two dimensional systems have been 
relatively more abundant, the associated rates of various heat and mass transfer processes 
(melting, dissolution etc.) in such systems have not been adequately, theoretically investigated. 
A typical problem in this arena is the melting of solid additions in two dimensional, two phase, 
turbulent flow systems (axisymmetrical gas stirred ladles). Thus, despite melting in 
axisymmetrical gas bubble driven systems been a relatively well studied phenomena, from 
experimental point of view, not much mathematical modelling has been earned out 
[10,11,12]. In particular, the usefulness of mathematical models to investigate melting 
phenomena in eixisymmetrical systems remains yet to be confirmed, particularly with 
reference to experimental data in high temperature systems. 

In this thesis a 2-D, steady state, two phase, turbulent flow model was used to predict 
the distribution of flow variables and turbulence parameters in axisymmetric gas stirred vessels 
for both aqueous and high temperature systems. The hydrodynamic predictions were 
subsequently, coupled with a previously reported heat transfer correlation [13] to predict 
melting rates of spherical shaped additions in the gas stirred systems. The predictions from the 
.model were compared with the experimental observations reported in literature on both aqueous 
,ind high temperature systems [14,15]. Alternative to this, a simple, yet effective, macroscopic 
)rocedure based on the ' Plume Model' [16] has been suggested to predict the melting rates of 
olids in gas stirred systems. Furthermore, two different but popularly used heat transfer 
Direlations were tested to predict melting rates in order to assess their relative merits and 
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demerits. 


In the second phase of the present work, the two dimensional code was modified for 
modelling a 3-D, steady state, two phase, turbulent flow situation (incorporating the two 
equation A:- 5 model of turbulence ) in cartesian co-ordinate system. The three dimensional code 
developed was extensively tested against available benchmark solutions. Subsequent to these, 
appropriate boundary conditions were incorporated and the 3-D calculation procedure was 
applied to predict melt flow in a simple rectangular shaped tundish with no flow control 
devices [17]. 


1.3 Layout of the Thesis 

The main body of the thesis consists of six chapters. Chapter-2 deals with modelling 
of melting phenomena in axisymmetric gas bubble driven systems. In this chapter, the broad 
outline of the steady, two phase turbulent model that was used as a starting point for the present 
investigation has also been elaborated upon. In Chapter-h the development of a three 
dimensional, steady state, two phase turbulence model starting from its 2-D counterpart has 
been presented. This included the governing equations, salient features of the formulation and 
model validation against benchmark solutions. Finally, in Chapter-4, the modelling of melt 
flow in metallurgical tundish systems are presented. Chapter-5 summarises the general 
conclusions derived from the present work and the thesis concludes with the recommendation 
for future work in Chapter-6. 
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CHAPTER -2 


MODELLING OF MELTING RATES IN AXISYMMETRICAL GAS 

STIRRED SYSTEMS 


As mentioned in the introduction ( Chapter-l ), an existing steady state, two phase, two 
dimensional turbulent flow model [ 18 ] was used as a starting point towards developing a three 
dimensional calculation procedure. In this chapter, the two dimensional turbulent flow model 
[ 18 ] is applied to investigate melting of solids in aqueous and high temperature axisymmetrical 
gas stirred systems. Once prediction from the existing CFD package is validated against 
experimental observations, reliability of the procedure [ 18 ] can be ascertained and hence, 
modification can be attempted with confidence. 


2.1 Introduction 

Melting is integral to many pyro-metallurgical operations carried out in both ferrous and 
non-ferrous industries. Large tonnage of scrap, ferro-alloys and de-oxidiser elements are added 
during the production of common metals and alloys. Melting of these additions consumes 
appreciable energy and hence, affects the process economics significantly. Consequently, there 
is considerable interest in identifying and understanding the phenomena at work during melting 
and subsequent dissolution of solid additives. 

Melting is the result of heat transfer by free and / or forced convection from the bulk 
liquid to the surface of the solid. The flux of heat to the solid surface is expressed eis : 

q-hiT,-T_) ( 2 , 1 ) 
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The corresponding rate of melting (say for a spherical shaped addition) can be readily derived 
from an appropriate heat balance and represented according to : 


dt 


h 




( 2 . 2 ) 


Since and are typically constant for a given system, consequently, as seen from Eqs(2. 1 ) 
and (2.2), the rate of melting, apart from the thermophysical properties, depends solely on the 
value of heat transfer coefficient, h. The heat transfer coefficient, intum, depends on the flow 
conditions, thermo-physical properties of the bulk phase etc., and must be known a-priori so 
as to estimate the melting rates effectively via Eq.il.l). Consequently, detailed information on 
fluid flow is required for the analysis and interpretation of melting phenomena. 

Since rates of various transport processes can be greatly influenced by manoeuvring the 
flow conditions within a reactor vessel, consequently, submerged gas injection today has 
become dereguer to metal refining operations carried out in ladles, furnaces and transfer 
vessels. The gas rising as the two phase plume to the free surface induces a turbulent 
recirculatory motion of fluid and thereby, as is well known, exacerbates melting and / or 
dissolution of solid additives. In the subsequent section, the detailed calculation [18] procedure 
is first outlined and it is subsequently demonstrated how melting rates of solid additives can 
be predicted from the first principles. 


2.2 Previous Work on Melting of Solids in Gas Stirred Baths 

The subject of melting and dissolution of solid additions in gas stirred vessels 
has received considerable attention over the years. Many studies have been earned out in 
stagnant melts while relatively few studies have been made in stirred vessels. Thus, numerous 
studies on the melting of solids in metallurgical melts have been reported in the literature [10- 
12,14,15,19,20]. Here, however, studies pertaining to gas stirred bath have been considered and 
reviewed very briefly. 
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In general, two different approaches have been adopted by the researchers to 
theoretically deduce the relevant flow velocities and hence, the melting rates (e.g., Nu = 
f{Re,Pr)). These for example included, 

i. a simple macroscopic energy balance [10,14] and, 

ii. the numerical solution of the turbulent Navier-Stokes equations [14,1 1,12,20]. 

Szekely and co-workers [10] were among the first to investigate melting of ice cylinders 
in an aqueous model of gas stirred ladle system. On the basis of numerical solution of the 
turbulent Navier-Stokes equations and experimental observations, the following correlation was 
proposed: 

Nu = 0.iiReTuf^(Pry^ (2.3) 

Subsequently, Taniguchi and co-workers [14,15] carried out melting experiments in the 
plume region of aqueous and high temperature gas stirred baths. These authors [14,15] 
demonstrated that a modified form of Ranz-marshall correlation viz., 

Nu-2-. (0.4 * 0.06 ^ (2.4) 

predicts melting rates fairly realistically and provides estimates that agree well with the 
experimental measurements. To this end however, Taniguchi et.al. [14,15] adopted a simple 
energy balance approach and proposed an expression, for estimating the characteristic velocity 
of the system (viz. the rise velocity of the two phase plume). 

-^(nr^p^v)v^ “ ITg + 0.06 (2.5) 

On the basis of such , Reynolds numbers , Re^ , {Eq.(lAy) and subsequently the Nusselt 
numbers were estimated. Predicted results were validated using experimental observations 
derived from high temperature and aqueous systems ( see Fig.(2.\)) None the less, as 
mentioned in detail towards the end of this chapter, the formulation of Taniguchi et.al. [14,15] 
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suffers from serious drawbacks. Consequently, an attempt has been made to simulate the 
experimental observations of Taniguchi et.al. [14,15] from a fundamental standpoint with the 
aid of a steady, two phase, two dimensional calculation procedure [18]. 


2.3 Mathematical Modelling of Fluid Flow in Axi-Symmetric Gas Injection Systems 

As mentioned already, melting rates in turbulent recirculatory flow system can be 
estimated provided the detailed information on flow fields in the system is known a-priori. 
Consequently, it is worthwhile to examine the flow calculation procedure [18] in some detail 
before attempts are made to predict melting rates. 

There are two different computational procedures typically been adopted to 
mathematically describe the hydrodynamics of two phase buoyancy driven systems and these 
include: 

i. models derived on the Eulerian frame of reference which are further classified into two 
groups viz., 

a. the quasi single phase models and, 

b. the two phase models 

ii. models derived on the Langrangian-Eulerian frame of reference. 

Traditionally, the relatively simple quasi-single phase calculation procedure has been 
more popular. However, two phase calculation procedures are relatively more rigorous and 
generalised. In the present study, as mentioned already a steady state, two phase, combined 
Langrangian-Eulerian approach has been adopted to model the flow phenomena in the 
axisymmetrical gas stirred system. As against a Eulerian two phase calculation procedure, the 
present method is computationally less demanding [18] . In this, the motion of the gas phase 
is computed numerically in a Langrangian frame in the form of steady stream of bubbles, while 
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the liquid motion together with the associated level of fluid turbulence were calculated by the 
routinely applied Eulerian scheme. 


2.3.1 Modelling assumptions 

To mathematically model submerged gas injection operations in axisymmetric vessels 

using a combined Langrangian-Eulerian approach, following idealisations and approximations 

were made: 

i. The presence of any overlying second phase was ignored and the melt-air surface was 
assumed to be flat. 

ii. The melt was assumed to be isothermal. 

iii. Following axial symmetry in the present geometry of gas injection, variations of flow 
properties along the angular direction were ignored. 

iv. Discrete mono-size bubbles were assumed to be form at the lance tuyere tip. Thus the 
size of the bubbles formed at the tip was known a-priori. From the principle of volume 
continuity the corresponding bubble frequency was calculated. 

V. Bubble-bubble interactions were assumed to be negligible and hence, drag coefficient 
correlation, applicable to the single bubble-fluid system was applied for the appropriate 
range of Eotvos and Reynolds numbers. 

vi. Bubble within the upwelling plume contribute to the generation of turbulence and this 
intum affects the bubble motion. 

vii. The flow is essentially Newtonian and incompressible and finally, 

viii. The flow is steady, turbulent and recirculatory. 


2.4 Description of the Steady State, Two Phase, Two Dimensional Turbulent Flow Model 

A brief outline of the combined Eulerian-Langrangian two phase, two dimensional 
turbulent flow model developed earlier by Mazumdar and Guthne [18] is presented below. The 
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governing equations together with their associated boundary conditions applicable to 
axisymmetrical gas stirred systems are summarised below in terms of the polar coordinate 
system. 

2.4.1 The liquid phase flow equations 

The governing equations in terms of the cylindrical polar co-ordinate system ( r, 6, z) 
chosen, can be represented as: 

Equation of continuity: 

^ (p, ^ (Pi a, ^ v) = 0 (2.6) 


Equation of motion in the axial direction : 


- -^(P;“/ 'w) = -a,-^ * 
az r dr dz 


^ d . du. 13/ 

2— (a,P — ) - (rn,a,— ) 

dz dz r dr dr 


Id, dv. E, 

(rii a, — ) + F 

rdr^ ‘dz’ ‘ 


(2.7) 


Equation of motion in the radial direction: 


d 13 

— (p,a, uv) * — — (p^a, rw) = -a 
dz r dr 


dp d . 3v. 2d. 3v. d . du. 

dz dz dz r dr dr dz dr 


2vp,a, 


r 


( 2 . 8 ) 


2.4.2 The liquid phase turbulence model 

In the present study, a modified form of standard coefficient k-e model of Launder and 
Spalding [21] was used. The governing equations for turbulence kinetic energy , k and its 
dissipation rate , e , can be represented in the cylindrical polar co-ordinates ( r, 6, z) as: 
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Equation of turbulence kinetic energy, k : 


02 r dr r dr dr dz dr 


(2.9) 


Equation of dissipation rate of turbulent kinetic energy, e : 


C,G*e 


y(p,a^e ) . -|-(p,a/Te) = lA(ra,^^) . ^(a ^ «/(^^ ' ^zPiT ^ At) (2. 10) 
dz r dr r dr a dr dz a dr k k k ^ ’ 


In Eqs(2.9) and (2.10), G^, is the volumetric rate of turbulence production, and is defined as: 


^ ,dv,^ ,v.\ ,du 

G* = Hj. [2[(— ) * (— ) *(-)]*(—*—)] 
dz dr r dr dz 


( 2 . 11 ) 


The eddy viscosity , , appearing in Eqs(2.7) and (2.8) is defined as : p^= pj-+ jUi ,m which. 


Pr = 


C,P,k^ 


( 2 . 12 ) 


The five empirical constants of the k-e model viz. ,Ci and were assigned to 

their standard values [21] . These are summarised in Table-2.1. 


Table-2.1 . Values of constant used in a k-e model. 


C, 

^2 





1.43 

1.92 

1.00 

1.30 

0.09 


2.4.3 Equation of motion of the gas phase 


As mentioned already, the motion of the gas phase has been carried out via a 
Langrangian approach, in which the trajectories of the individual bubbles are stochastically 
determined in space by solving two differential equations in time. The appropriate form of 
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Newton's second law of motion that describes the bubble motion, can be represented along the 
two co-ordinate axis , r and z, according to: 


Axial direction, z : 


(1 . 0.5-^) 




dt 


- U) . . V^] . (1 - ^)g 

P, az ar 


(2.13) 


Horizontal / Radial direction, r: 


(1 


0 . 5 —)^^ 
P, dt 


-^C^e - V) 


P, az dr 


(2.14) 


The two corresponding kinematic relationships, which define the trajectories of bubbles within 
the two dimensional flow domain are: 


dz 

It 


B,x 


(2.15) 


and , 



(2.16) 


The instantaneous drag coefficient, C^,, in the Eqs(2.\3) and (2.14) was evaluated using an 
experimentally determined correlation [22] for oblate spheroid / spherical capped bubbles: 


2:^22 . < Re < 5000 ) 

[— * 0.235] 


(2.17) 


Furthermore, Re in £'^j'(2.13) and (2.14) is the Reynolds number, and is estimated on the basis 
of absolute magnitude of the relative velocity vector, , i.e. 
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(2.18) 


U,-U 
Re=p,dJ-^ 

Similarly in Eqs(2.\3) and (2.14), U and V are the instantaneous components of fluid motion, 
and were estimated by summing up the relevant mean velocity components ( u and v) with the 
corresponding fluctuating velocity components assuming turbulent fluctuations are essentially 
isotropic and obey a Gaussian probability distribution function. On the basis of such, 

i7=„,T(y)^ (2.19) 

where t|; is a normally distributed random variable and varies between -1 and 1 . Thus during 
the rise of the bubble through the liquid, the velocity fluctuations were estimated stochastically 
and were assumed to remain constant for a time period, equivalent to the eddy life time, r . The 
life time of an eddy, during which a characteristic velocity fluctuation was assumed to persist, 
is estimated from the turbulence theory according to [23]; 

X = 0.195 (2.20) 


2.4.4 Estimation of gas voidages 

Figure (2.2) illustrates a typical section of the grid network and a bubble trajectory 
intercepting the Eulerian grid system. For a scalar control volume the corresponding gas 
voidages can be readily estimated from the principle of volume continuity, according to : 

volime of gas s 

cc i 2m \ I 

‘ volume of control volime 2A ^ 

These are summarised elsewhere in detail [18] and consequently, not reproduced here. Once 
Ug values have been determined for each of the control volumes, corresponding liquid volume 
fraction, a , , can be easily determined, from + a; =1.0. 
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Fig.{22) 



« scalar nodes o bubbles 


: Schematic of a section of Eulerian grid network intercepted by a bubble 
trajectory. 
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2.4.5 Estimation of turbulence production ( P* ) due to bubbles and interphase friction 
forces i F. ) and {F^ ) 


Bubbles in a two phase region contribute significantly to the total turbulence generation 
[24] . The additional turbulence produced by the bubbles can be considered to be equivalent to 
the shear work performed on the liquid by the bubbles and therefore, the corresponding 
volumetric flow rate of turbulence generation, , can be represented as : 




Q. 


c 








( 2 . 22 ) 


where, Vg is the resultant bubble velocity and Uf is the resultant fluid velocity. Q in Eq.(2.22) 
is a constant and takes a value between 0 and 1. The value of Q depends on the physical 
characteristic of the two phase plume ( e.g., bubble population, bubble size etc.). The constant 
Q is normally determined by trial and error. 


The momentum exchange terms F. and F^ can be evaluated assuming that drag forces 
experienced by bubbles act with equal magnitude but opposite directions in the liquid phase. 
The expressions for F. and P). can be represented in the compact tensorial notation as : 


F = 


Q. 


NAV. 




4d: 


( 2 . 22 ) 


2.4.6 The boundary conditions 

The boundary conditions applied to the numerical solution of the set of governing 
equations for fluid flow and turbulence in axisymmetrical gas stirred systems are represented 
schematically in Fig. (2.3) and are as follows; 

At the axis of symmetry {r = 0 , 0 <z< L) 
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.. 0 ; f = 0 ;|. 0;|.0 

dz oz cte 


X. 


o 

II 


^1 




o 

o 

11 

> 


o 

II 




ro 


O 

II 

Ci) 

6" 

II 

o 

II 

> 

o' 

II 


• I I I I I I I I I I I I I , ■ 

Z = Q J 


wall functions applied 
on these nodes 


Gas injection 
nozzle/ plug 


Fig.(23) : Schematic of the vessel geometry and the corresponding boundary 

conditions applied towards numerical simulation of axi-symmetric gas 
stirred ladle system. 
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At the side walls (r = R,0<z<L) and bottom surface {z - 0 ,0 < r < R) 

u^v=k=e=0 

Close to the rigid boundaries ( viz., the walls and the bottom surface ) where the variations in 
the flow properties are steep, the momentum ( u and v ) and the scalar field ( k and e ) have 
consequently been modelled using the routine wall function procedure [21,32] i.e. logarithmic 
law of wall for the velocity component parallel to it (see Chapter-4 for details). In addition to 
these the appropriate gas inlet velocity was applied as the required initial condition through 
£:^s.(2.13) and (2.14). 

2.4.7 The numerical solution procedure 

The basic structure of the present combined Langrangian-Eulerian procedure has been: 

i. The total flow is equally divided between a preconsidered number of trajectories. On 
the basis of such, flow rate (Q,) per trajectory is estimated. 

ii. The equations for bubble motion are solved for each of the trajectories and then 
superimposed on the Eulerian grid. From these, overall distributions of a/ , F. and F,. 
etc are estimated. 

iii. Estimated <ar, , F. and F) are incorporated into the equations of the continuous phase 
( viz., the liquid ) and solved iteratively by embodying the SIMPLE procedure of 
Patankar and Spalding [26]. 

iv. At the end of a predetermined number of iterations of the liquid phase equations, the 
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trajectories of the bubbles are recalculated using the prevalent flow and turbulence 
fields. Based on this the distribution of a, , F. and etc. were updated. The above 
steps are repeated till the converged solution is not obtained. 

The flow circuit of the computer program applied to the present investigation has been 
illustrated in Fig.(2A). 


2.5 The Heat Transfer Model 


Once the flow velocities and turbulence parameters are predicted as a function of 
operating conditions, one of the many empirical expressions available in literature can be 
applied to quantify heat transfer between a submerged solid sphere and a continuous phase. For 
the system configuration considered in this work ( see Fig.Q ..\) ) , the relationship proposed 
by Whitaker [13] was adopted. 


Nu-2 = ( 0.4 




0.25 


(2.24) 


In £'^.(2.24), Re ^ , is the object Reynolds number that embodies the initial diameter of the solid 
sphere and the characteristic velocity prevalent in the neighbourhood of the solid sphere. It is 
instructive to note here that £^.(2.24) was originally developed for heat transfer between a 
heated metal sphere and a bulk single phase fluid ( air ). As discussed elsewhere [13-15] , the 
above correlation was found to be applicable to turbulent, two phase flow systems as well. 


2.6 Results and Discussion 

To assess the predictive capability of the model developed in the present study ,the 
experimental study reported earlier by Taniguchi et.al. [14,15] was applied as the reference. 
Thus, for the given vessel geometry (see £/g.(2.1)) and gas flow rates, the flow distributions 
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Fig.{2A) : Flow diagram of the two dimensional, two phase turbulent flow model 
[18]. 
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were estimated first through the numerical solution of Eqs.(2.6) through (2.23). Subsequent 
to this the Re)mold numbers, Re^ , were calculated ( Table-2.2) via. the following expression, 
embodying the predicted velocity prevalent in the vicinity of the solid spheres viz.. 




D p, 


1 ^/ 


(2.25) 


Incorporating such estimates and the relevant thermo-physical properties ( Table-2.2) in 
£^.(2.24), Nusselt numbers were estimated. In Fig.(2.5) a log-log plot between the 
experimental Nusselt number [14,15] and the present estimates of the Reynolds number (see 
Table-2.2) is shown together with the theoretical line for £^.(2.24). There the dependency of 
the experimental Nusselt number on the numerically predicted Reynolds number appears to 
have been simulated reasonably well by the present mathematical model. 


Table-2.2: Characteristic of the experimental work of Taniguchi et.al. and the present 
estimates of the liquid rise velocity and Re^ in the aqueous and high temperature systems : 


(a) Aqueous Model 


Experimental set 

up 

Temp, 

7;,(K) 

Fx 10®, 

(Nm^/s) 

0x10®, 

(m’/s) 

Mean Liquid rise velocity 

in the vicinity of the sphere, 

{m/s) 

Estimated Re^ 

Numerical 

Prediction 

Eq.(2.28) 

Numerical 

Macro- 

scopic 

L = 0.36 m 








R = 0.075 m 

293 

1.07 

1.03 

0.063 

0.082 

2259 

2926 

D = 0.036 m 








;/*= 1.02x10’ 

293 

7.38 

7.11 

0.155 

0.157 

5558 

5579 

//o=1-83x10-’ 








Pr = l.\5% 

293 

29.51 

28.41 

0.301 

0.247 

10758 

8853 


@ all units in S.I.system 
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Fig.(2.5) : A log-log plot between the experimental Nusselt number [ 1 4, 1 5] , 
(Nu-2)/{Pr°"* (ni, /juof'^^ } the present estimates of the Reynolds 
number, Rcp ,vis a vis the line for £'^.(2.24). (Deduced through numerical 
simulation ). 
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(a) High temperature Model 


Experimental set 

up 

Temp 

,T„ 

(K) 

Fx 10®, 
{Nm^/s) 

QxlO®, 
{m^/s ) 

Mean Liquid 

in the vicinity 

(jTV 

Numerical 

Prediction 

rise velocity 

of the sphere 

's) 

Eq.(2.28) 

Estimat( 

Numerical 

sd Rbq 

Macro- 

scopic 

L =0.2 m 

953 

0.98 

3.01 

0.067 

0.124 

2206 

4083 

R = 0.04 m 

953 

3.08 

9.47 

0.132 

0.182 

4347 

5993 

D =0.024 m 

953 

8.83 

27.13 

0.196 

0.258 

6454 

8496 

//553= 1.72 X 10-^ 

973 

0.98 

3.08 

0.069 

0.125 

2383 

4318 

//973= 1.64 X 10"^ 

973 

1.58 

4,97 

0.092 

0.147 

3178 

5077 

/’r ,53 = 0.0211 

973 

2.08 

6.54 

0.112 

0.161 

3869 

5561 

P/-,53 = 0.0197 

973 

3.08 

9.68 

0.142 

0.183 

4912 

6321 


@ all units in S.I.system 


Numerical results as are well known, depend on various assumptions and 
approximations applied to the mathematical model. In the present study the predicted results 
are expected to be sensitive to the precise choice of values of the bubble diameter dj,, drag 
coefficient Cp, the empirical constant Q , as well as to various numerical parameters such as 
grid size, convergence criteria etc. To access the affect of these parameters on the computed 
results, extensive numerical experimentation was carried out for the aqueous system at the gas 
flow rate of 29.51 x lO"^ Nm^/s. To this end several grid systems viz., 14 x 8, 18 x 10, 24 x 16 
were tested to arrive at the grid independent solution. The liquid rise velocity in the vicinity of 
the solid sphere obtained by the three grid systems are shown in Table-23. 

These clearly indicate that refining the grid beyond (18 x 10) did not bring in any 
significant changes in the -predicted results. Similar trials with different convergence criteria 
indicated that the absolute value of srrai on residuals on u, v, p, k and e over the whole flow 
domain has to fall below 5x10'^ (S.I.units) in order to arrive at a converged solution. All the 
results provided in this text were derived on the basis of a 18 x 10 grid system and the above 
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mentioned convergence criterion ( for definition see Chapter-A). 


Tcible-2 .3 : Sensitivity of the predicted axial velocity to the various grid configurations applied. 
( Aqueous system, L = 0.36 m,R = 0.075 m,Q= 29.51 x lO"* NmVs ) 


Grid System, 

(zxr) 

Mean rise vel, in the 

vicinity of sphere, 

(m/s) 

Estimated 

Estimated 

Nu 

Estimated 

log 

14x8 

0.23 

8616 

120.4 

2.08 

18x 10 

0.30 

10758 

136.5 

2.14 

24x 16 

0.29 

10411 

134.0 

2.13 


Considerable uncertainties can be associated with the value of the bubble diameter 
which forms an integral input to the present mathematical model. In a gas stirred system a wide 
spectrum of bubble size can be envisaged in the plume. Consequently, a fair degree of 
idealisation has been involved while applying a single characteristic value of bubble diameter 
for such systems. Hence, the bubble diameter value was varied arbitrarily in the computational 
procedure within a range of ± 30 % from the reported experimental value [14,15]. While 
varying the bubble diameter value for a given gas flow rate, due care was taken so that volume 
continuity remains satisfied. Numerically predicted axial velocity component as a function of 
three different bubble diameters is shown in Fig. (1.6) for a dimensionless height , z/H= 0.5. 
There, as seen, the axial velocity within the plume region appears to be sensitive to the value 
of the bubble diameter, while outside the two phase region, predicted axial velocity seems to 
remain practically unaffected by such variations. For a given flow rate, finer sized bubbles 
increase the surface area to volume ratio and therefore lead to a higher drag force ( oc 1/d,^ 
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2/H =0.5 


- 30 % 

db (=0.010 m) 

+ 30 % 


\ 

, \ 

V \ 

\ \ 

• \ 

■•.Xn 


Dimensionless Radial Distance 

Fig.{2.6) ; Senstivity of the predicted vertical velocity component to the diameter of 
monosize bubbles (d^) assumed to be forming at the lance tuyure tip. 




) which in turn leads to higher fluid rise velocity, particularly in the plume region. 


Similarly, the drag coefficient (C^,) appearing in £’^5(2.13) and (2.14) was also varied 
in the range of ± 10 %. It is important to mention here, that drag coefficient is a function of the 
bubble diameter. Consequently, uncertainty in the value of the bubble diameter is also likely 
to introduce some uncertainty ( of course, to a smaller extent ) in the value of the estimated drag 
coefficient. As illustrated in Fig. (2.7) the sensitivity of the predicted axial velocity component 
to different drag coefficient values was found to be minimal in the bulk ( outside the plume 
region ) of the liquid. However, in the plume region, the rise velocity was found to vary 
appreciably, almost to a factor by 2. Since a higher drag coefficient leads to a higher drag force, 
consequently, one can anticipate that the rise velocity of the liquid in the two phase plume 
region will increase with any increase in 'C^' value. 


Table-2.4: Variations in the estimates of Nusselt number as a function of the bubble diameter 
(dj, drag coefficient (C^^) and coefficient of turbulence generation { Q ). 


Parameter 

Numerical Value of 

Mean rise vel. in the 

Estimated Nu 


the Parameter 

vicinity of sphere, (m/s) 

number 

Bubble Diameter 

0.007 m 

0.36 

159.5 

(d,) 

0.010 m 

0.30 

136.5 


0.013 m 

0.19 

105.8 

Drag Coefficient 

2.40 

0.26 

126.1 

(Co) 

2.67 

0.30 

136.5 


2.94 

0.35 

148.8 

Empirical constant. 

0.0 

0.30 

136.5 

Q 

0.5 

0.26 

126.1 


0.7 

0.22 

114.8 


0.9 

0.09 

105.8 
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The coefficient 'Q' appearing in Eq.(2.22), as mentioned already, is a function of the 
bubble population, shape and size etc. and hence, is difficult to precisely quantify in such gas 
stirred system. To ascertain the influence of this variable on the predicted results, four 
different but plausible values of 'Q' ( 0, 0.5, 0.7 and 0.9 ) were applied to the calculation 
procedure. Results thus obtained have been illustrated in Fig. (2. 8). There, the axial velocity is 
seen to be dependent of some extent on the value of 'Q' particularly in the plume region. A 
similar trend was found for the number of bubble trajectories. Increase in the number of bubble 
trajectories was found to lead to a higher liquid rise velocity within the plume. It is instructive 
to note here that marked variations in the velocity within the plume region do not translate by 
the same order in the bulk, since the plume barely occupies 1 to 2 % of the reactor 
volume. 

The results presented so far appear to indicate that assumption on the bubble size, drag 
coefficient etc. embodied in the mathematical model are very likely to effect the predicted 
results, and hence the inferences drawn therefrom. The extent of such influence as shown on 
Figs(2.6) through (2.8), while is expected to be only marginal in the bulk phase region, is likely 
to be more pronounced in the two phase region of the vessel. The present work therefore, 
suggests that sufficiently accurate information on the bubble diameter f ri* j, the drag coefficient 
( Cq), and coefficient for turbulence generation for bubbles (Cf,) is needed, if fluid flow and 
the associated heat and mass transfer phenomena are to be predicted realistically via. such two 
phase calculation procedure. Fortunately, however as the melting rates are a weak function of 
velocity, Nu « u''^ (see £'^.(2.24)), such uncertainties in the predicted rise velocity are likely to 
have relatively small influence on the predicted melting rates of solids as illustrated in Table- 
2.4. 

2.7 An Alternative Macroscopic Formulation 

In the absence of any computational facility, it is worthwhile to examine the adequacy 
of a simple alternative procedure for predicting melting rates in such systems. In their extensive 
study of the melting phenomena Taniguchi at.al. [14,15] proposed a simple energy balance 
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Axial Velocity, m/s 



Fig.(2.S) : Senstivity of the predicted verticle velocity component to the value of 
coefficient of turbulence generation by bubbles (Q). 
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relationship for estimating the liquid rise velocity in the two phase gas-liquid plume region, 
e.g.. 


^(■Kr^p,vW-W,^OMW^ (2.26) 

In Eq. (2.26), Wg and are respectively, the rate of buoyancy and the kinetic energy input to 
the system and were assumed known from the operating conditions viz., vessel radius, gas flow 
rate, nozzle diameter etc. Similarly v is the mean upward flow in the two phase region (see 
F/g.(2.1)) and r is the average radius of the bubble plume. The value of r was determined 
experimentally by Taniguchi et.al. and embodied in Eq.(2.26) for estimating the mean upward 
velocity , v , in the two phase plume region. 

Numerous experimental studies reported to date on the gas stirred ladle system appear 
to indicate that the measurement of the gas-liquid plume parameters is somewhat uncertain 
owing to the dynamic nature of the bubble plume. Such two phase plume, in general exhibits 
long and short term wondering [27] and thus makes it difficult to measure r realistically. In 
industrial gas stirred reactors, reliable measurements are even more difficult owing to the 
opacity of melts. Assuming that the kinetic energy input is negligible Eq. (2.26) leads to : 

Kr^v^^p,gQL (2.21. a) 

where, the rate of buoyancy input Wb = p, gQL. Hence, for a given set of operating conditions 
Eq. (2.26) can be conveniently written as: 

rccv** (2.27) 

Consequently, any uncertainty in the measured plume radius ' r ' is likely to effect the 
estimated flow velocity ' v ' more pronouncedly. As a consequence of these, an alternative, 
physically more plausible formulation is called for the estimation of the mean upward flow and 
the associated rate of melting in the two phase plume region of the gas stirred systems. 

Since the specific potential energy input rate of the experimental configuration of 
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Taniguchi et.al. is similar to those considered by Mazumdar and co-workers [28], consequently, 
the ' plume model ' derived earlier [28] for ladle metallurgy steelmaking operations can be 
conveniently extrapolated to their experimental study. Equating the potential energy input rate 
to a fraction of the rate of turbulence energy dissipation in the gas stirred systems, Mazumdar 
et.al. [28] demonstrated that the rise velocity of the up-welling gas-liquid plume mixture can 
be estimated in terms of the operating variables L, R and Q according to, 

(2.28) 

in which, L is the depth of the liquid, R is the vessel radius and Q is the gas flow rate 
corrected to the mean height and the temperature of the liquid. Details of the derivation can be 
found elsewhere [21], and is therefore not reproduced here. The adequacy of the ' Plume Model' 
has already been demonstrated independently by several researchers [29,30,31]. Consequently, 
it is proposed that the mean upward velocity in the two phase plume region be estimated via 
£^.(2.28) in lieu of Eq.(2.26), and applied in conjunction with an appropriate heat transfer 
correlation for estimating the melting rates of solids in the two phase plume region. Thus Up, 
the plume rise velocity estimated a-priori via £'q.(2.24) in terms of the key operating parameters 
(viz., L, R and Q) can now be embodied in an appropriate heat transfer correlation to deduce 
the heat transfer coefficient (viz., Eg'. (2.24)) in the two phase plume region of the gas stirred 
system. 


Thus considering, Re^, = (p,UpD)/ ///predictions on melting rates were made for the 
experimental observations of Taniguchi et.al. [14,15] were quantified on the basis of £^.(2.24). 
As mentioned already, the experimental configurations and the values of relevant 
thermophysical properties together with the present estimates of plume velocities as a function 
of the gas flow rates applied have been summarised in Table-2.2. In Fig.(2.9) a log-log plot 
between the experimental Nusselt number [14,15], (Nu-2)/{Pr°"' (pi, / A ^nd the 

corresponding estimates of the Reynolds number, Re^, (pjUp D )/ p,)is shown together 
with the theoretical line for £'^.(2.24). There very reasonable agreement between estimates and 
the experimental measurements is readily apparent. 
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Fig.(2.9) : A log-log plot between the experimental Nusselt number [14,15], 
(Nu-2)/{Pr°'* (Hk } and the present estimates of the Reynolds 

number, , vis a vis the line for Eq.(2.24). ( deduced through the 
macroscopic modelling Eq.(2.2S). 
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Alternatively to £'^.(2.24), the following heat transfer correlation has also been applied 
to assess the adequacy of the macroscopic formulation further e.g., 

A^m - 2 = [ 0.06 Pr ° ] (2.29) 

It is important to note here that in an earlier study Iguchi et.al.[20] have demonstrated 
that Eq. (2. 29) holds good for gas stirred system and produces estimates of Nusselt number that 
were in reasonable agreement with their experimental observations as well as the ones reported 
in literature. In the plume region of a gas stirred system, the intensity of turbulence is known 
to be fairly large [24] in comparison to those in the bulk phase region. To this end, the 
experimental work of Sheng and Irons [24] suggests that wit hi n the plume of an aqueous gas 
stirred system, the intensity of turbulence can be ascribed with reasonable certainty to a 
characteristic value of around 0.5. Based in such information Eq. (2.29) can be simplified for 
the two phase region as: 

(2.30) 

On the basis of £^.(2.30)the Nusselt number and / or the heat transfer coefficients were 
re-estimated in terms of the operating variables. In £/g.(2.10), a log-log plot between the 
experimental Nusselt number [14,15], (Nu-2)/Pr°^^ and the estimated Reynolds number, 
(=(Pi Up D )/ Pi), is shown together vrith the theoretical line for Eq. (2.30). There, reasonable 
agreement between experimental measurements and present estimates is at once apparent. 
Results presented in £/g.(2.10) evidently indicate that the present formulation is indeed as 
effective as the differential model predictions [20] for estimating the melting rates in the two 
phase plum region of the gas stirred baths. 

As a final point the present study also indicates that the Eqs (2.24) and (2.29) are equally 
effective for predicting of melting rates in the gzis stirred bath. However, Eq. (2.29) embodies 
an additional dimensionless group i.e., the intensity of turbulence, Tu , in addition to the object 
Reynolds number, Re^ The distribution of this parameter should be known a-priori so as to 
estimate the heat transfer coefficient via. Eq. (2.29) particularly, when the addition is melting 
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Fig. (2.10): A log-log plot between the experimental Nusselt number [14,15], 
(Nu-2)/ Pr° and the present estimates of the Reynolds number, Re ^ , 
vis a vis the line for £'^.(2.30).( deduced through the macroscopic 
modelling Eq.{2.2S). 
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outside the plume in the bulk liquid. Analysis of melting under such conditions will therefore 
warrant detailed computation of the turbulent flow field through numerical calculation. 
Nevertheless, it is well known that while the flow variables ( w, v etc.) in such systems can be 
estimated with reasonable certainty, considerable uncertainty still remains as far as prediction 
of various turbulence parameters (p^ , k etc.) are concerned. Evidently, therefore some 
uncertainty is likely to be associated with the numerical estimates of Tu, the intensity of 
turbulence. In view of such limitations associated with the hydrodynamic modelling of gas 
stirred systems, it appears that £^.(2.24) is likely to provide somewhat superior and more 
reliable estimates of melting rates in gas bubble driven systems, than Eq. (2. 29). 
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CHAPTER -3 


DEVELOPMENT OF A STEADY STATE, THREE DIMENSIONAL 
TURBULENT FLOW CALCULATION PROCEDURE 


In the previous chapter, the adequacy of the two phase computational procedure 
developed by Mazumdar and Guthrie [18] was examined. It was demonstrated implicitly that 
their model produces estimates of flow and turbulence parameters in the gas stirred system 
that are in reasonable agreement with the experimental observations reported in the literature. 
In this chapter the salient features of upgradation of the CFD package [18] to steady state, 
three dimensional, two phase, configuration is presented. It needs to be mentioned here that 
the 3-D calculation procedure has been developed for cartesian coordinate system, while the 
calculation procedure of Mazumdar and Guthrie has the capability of carrying out 
computations in both cartesian and polar coordinate system. 


3.1 Introduction 

The degree of complexity required to mathematically model different metallurgical 
processes varies from one to another. While, operations like axisymmetric gas injection, 
axisymmetric pouring and drainage, can be effectively modelled using a two dimensional 
turbulent flow model, others, like melt flow in tundishes or in a R-H degasser require a three 
dimensional turbulent flow model for effective simulation. Needless to mention that transport 
processes in practically all metallurgical operations can be successfully modelled using 
transient, two phase, three dimensional models. To this end, several software packages are 
currently commercially available. Some of the popularly used software packages have been 
listed in Table-3. \. 
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Table-3.1: Popularly used software packages to study transport phenomena [2]. 


Name 

Application 

Technique 

PHOENICS 

fluid dynamics, Heat Transfer, Chemical Reaction 

finite volume 

TEACH -T 

fluid dynamics 

finite volume 

FLOW 3-D 

fluid dynamics 

finite difference 

FLUENT 

fluid dynamics, Heat Transfer 

finite difference 

NECTON 

fluid dynamics. Heat Transfer 

finite element 

NISA 

fluid dynamics. Heat Transfer 

finite difference 


However, the key problem lies in the definition of the boundary conditions and the 
ability to interpret the results in an intelligent manner. Hence, many researchers prefer custom 
built models from the first principles with consideration of elementary control volumes. 
Indeed, this later philosophy has been applied in the present investigation to develop the 
steady, two phase, three dimensional computational procedure. 

Thus, the present chapter deals with the development and validation of a steady state, 
three dimensional, two phase, turbulent flow model in cartesian coordinate system. The 
computer code developed has been subsequently used to study the melt flow behaviour in 
rectangular shaped steelmaking tundish system with no flow control devices ( Chapter-A). 


3.2 Governing Equations 

3.2.1 The liquid phase flow equations 

In terms of cartesian co-ordinate system {x,y,z ) , the governing flow equations ( for 
the liquid phase ) under turbulent flow conditions can be represented as: 
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Equation of continuity : 


-^(a pw) * —(a pv) ^ —(a pw) = 0 
dx dy dz 


(3.1) 


Equation of motion in x direction : 

d , . d , \ d , , dp d f du. d , du^ d , du. „ 

—(a, puu) - —(a, pwv) - —(a, ptw) = —(a, " ^(“/ 

dx 3y dz dx dx dx dy dy dz dz 

(3.2.a) 


in which, 


„ d , du, d , dv, d , .. Sw,. p 

“ '^'cbc dy ' ‘dx’ dz ' ‘dx’ ‘ 


(3.2.b) 


Equation of motion in the y direction : 

—(a pwv) » —(a, pw) + pvw) = - ^ 

dx ' dy dz dy dx dx dy dy dz dz 

(3.3.a) 


in which. 


„ d y 5Mn d / dV\ d y dw^ j-, 

■ s‘“' '*>> • P"' ■ p°‘ ' ’ 


(3.3.b) 


Equation of motion in the z direction : 


pwv) . ■|.(., pw) . |(«, pw) . - I • ^(a, P,.^) . M.^) . |(«, P.^) . s. 

(3.4.a) 


in which. 
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(3.4.b) 


C ^ 5 5 17 


Equations (3.2) through (3.4) are the well known Navier-Stokes equations. In these u, v, w 
and p (gauge with reference to local hydrostatic pressure) are respectively the time averaged 
velocity components and pressure. iXi is the volume fraction of liquid and F's are the various 
body forces acting on the fluid which may include buoyancy, drag etc. Finally, , appearing 
in Eqs.(7)2) through (3.4) is the effective viscosity and is derived from an appropriate 
turbulence model. 


The values of ai in Eqs (3.2) through (3.4) can be obtained from the solution of bubble 
trajectory equations already outlined in Chapter-2. Assuming, that there is no bubble motion 
in the z-direction (the lateral direction) the bubble trajectory equations can be used in its 
original two dimensional form, outlined in the previous chapter. Consequently, the drag 
forces generated by the bubble motion will exist only in the x and y directions i.e. F. will be 
equal to zero. Subtle modifications are warranted such as, the relative velocity vector between 
the bubbles and the liquid phase will now be required to be calculated from the three velocity 
components in the x, y and the z directions. Moreover, the two dimensional arrays of bubble 
motion need to be converted into the corresponding three dimensional arrays. Alternatively, 
an additional force balance equation can also be incorporated for the z coordinate axis and 
thus, components of body forces i.e., F^, F^, F.m. the Navier Stokes equations can be 
calculated by solving the three ordinary differential equations governing the motion of the 
bubble in the three dimensional rector. All calculations reported in this and subsequent chapter 
were derived assuming ai =1.0 (single phase calculation). 

3.2.2 The turbulence model 

In the present work the popular two equation, Model of Launder and Spalding 
[2 1 ] has been used for the estimation of turbulence parameters. The model incorporates an 


41 



equation for the conservation of turbulent kinetic energy, k{ - \l2(u'^ + v'^ +w'^) } and its 
dissipation rate, e { = -dk/dt } . The k and e transport equations in terms of cartesian co- 
ordinate system can be expressed as following: 


Conservation equation for the turbulence kinetic energy, k: 


^(P, ^(P, a-ivk) 

dx dy 


dz dx Oi^dx 


a, p-.dk. d, i^edk. 

dy Oj^dy dz " 


o^dz 


(3.5.a) 


in which. 


-5* = “; ( <5* - P^ ) 


(3.5.b) 


Conservation equation for the dissipation rate of turbulence kinetic energy, e : 


dx 


(p, a, ue) * ^(p, a, ve) + — (p, a, we) 


dy 


d, p.ae, d, P, ae. 
—(a, — — ) . —(a, — — ) 
dx a dx dy o ^ 


d , P. ae, 

—(a, — — ) 

dz o dz 


(3.6.a) 


in which. 


5 =«,[ 


C,eG^- 


Cpe^ 


(3.6.b) 


G* appearing in Eqs.Q.S) and (3.6) is the volumetric rate of turbulence generation and for a 
three dimensional cartesian co-ordinate system, can be expressed as: 


ax dy dz 


PrK 


du 

dy 


dvd ,dv 


dy 


r du 

'dz 




(3.7) 


The eddy viscosity , , appearing in Eqs(3.2) to (3.4) is defined as: pj-+ p/ , in. which. 


Pr 


C fi 


(3.8) 
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The five empirical constants of the k-e model viz., C, , C 2 , and were 

assigned to their standard values [21]. These have been already summarised in Table-2 . 1 (see 
Chapter-2). As seen from the above, flow and turbulence model equations are mutually 
coupled. 

3.2.3 The boundary conditions 

The boundary conditions are problem dependent and most of the metallurgical 
systems encountered include a free surface of liquid, symmetry axis ( if there is any) and the 
solid vessel walls . Boundary conditions for velocity components include no slip condition 
at the boundary walls (zero velocity) while across the free surface, zero shear is assumed to 
be transmitted. Gradients of all the velocity components at the axis of symmetry are taken as 
zero. Similarly, the values of k and e at the walls are set at zero. Zero gradient of k and e are 
applied at the symmetry plane and the free surface. Besides this, known values of the 
dependent variables are applied at the entry / exit locations depending on the problem. 

The near wall regions have two major characteristics that distinguish them from the 
bulk flow [21,32]: 

i. steep non-linear gradient in the flow properties, and 

ii. levels of Reynolds number are sufficiently low for local isotropy of small scale 
turbulence to prevail. 

Hence, the basic turbulence model becomes inapplicable in these regions and alternate 
representation must be sought. In the present study wall function procedure [21,32] has been 
used. The method is based on one-dimensional fluid flow, and assumes that the velocity 
profile in the near wall region is given by the universal logarithmic wall law which can be 
represented as : 
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(3.9) 


^-^\n{yE) 

U, K 


where u is the velocity parallel to the wall and u^{ = V (t^I p)) is the friction velocity, /ris 
the van Karmen constant and E is the roughness parameter typically taken as 9 for 
hydraulically smooth walls, is a dimensionless wall distance and is defined as : 


P) 

y = 

P 


(3.10) 


where y is the normal distance from the wall. This law has been applied to a point whose y' 
value is greater than 11.63. This is schematically shown in Fig. (3.1) . In the y" region 
specified above the Reynolds stresses are constant and the wall stresses, can be calculated 
as ; 

K-9C;^k (3.11) 


Further, in this region the convection-diffusion of (u,Uj ) is negligible. This implies that the 
rate of production of turbulent kinetic energy k is equal to its dissipation. Such considerations 
lead to: 


k 1 

(v/^) 


( 3 . 12 ) 


This relation is normally used as an boundary condition for k in the two equation k-e model. 
Similarly the boundary condition for e can be expressed as : 

6 = -^ (3.13) 

Ky- 

The preceding relationships are to be applied to a near wall node whose y^ is greater than 
11.63. 
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Fig.Q . 1) : Schematic representation of procedure for the near wall treatment of flow 
variable and turbulence parameter in the turbulent flow model. 
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3.3 The General Differential Equation 


The general structure of the relevant differential equations ( Eqs.{'^.\) to (3.6)) 
describing the conservation of mass and momentum and that of the k- e Model of turbulence 
indicate that all the dependent variables of interest obey a generalised conservation principle. 
If the general variable is donated by (p, the general differential equation can be written as: 

div ( p «(()) = div (rgrad^) i- (3-14) 


The three terms in the general differential equation, are the convection term, the diffusion 
term (/’is the diffusion coefficient ) and the source term, S^. The quantities /"and S are 
specific to a particular meaning of (p. The three terms in the general differential equation are 
the convection term, the diffusion term, and the source term. Depending on the general 
variable 0,'an appropriate meaning will have to be given to the diffusion coefficient T^and the 
source term 5. For the present set of governing equations this is summarised in Table-3.2. 

Table-3.2 : Definition of F, (p and S of £'^.(3.14). 


Conservation of 

0 

r 


mass 

1 

0 

0 

momentum 

u 

■HjH 

Eq.(3.2.a) 


B 


Eq.(3.3.a) 


Bi 


Eq.(3.4.a) 

turbulence kinetic 

B 


Eq.(3.5.a) 

energy 

B 



dissipation rate of k 

€ 


Eq.(3.6.a) 
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Hence, the procedure of casting any particular differential equation into the general 
form (£^.(3.14)) requires manipulation of the equation so that the convection and diffusion 
terms conform to the standard form. The coefficient of grad <p is taken as the expression for 
r, and the remaining terms on the right hand side are collectively defined as the source term, 
S. Since, all the equations for mass and momentum transport, and turbulence related 
phenomena can be thought of a particular case of general equation (£^.(3.14)) consequently, 
we need to concern ourselves with the numerical solution of £ 9 . (3. 14) only.. 


3.4 The Grid Arrangement 

A portion of the three dimensional grid is shown in £/^.(3.2). For the grid point P, 
points £ and W ( east and west respectively) are its x-direction neighbours, while N and S 
(north and south respectively) are the y-direction neighbours and T and B (top and bottom 
respectively) are the neighbours in the z-direction. The whole flow domain has been divided 
into a number of such non-overlapping control volumes. Staggered control volumes have been 
considered for the velocity components u, v and w. Fig.Q.S) represents the different staggered 
control volumes at which the discrete values of the variables are located. These values 
represent the averages over the respective control volumes. Scalar variables i.e. p, k,p^,p and 
e are stored at the grid nodes, while the velocities w, v and w are stored midway between the 
grid nodes, between the pressure which drives them. Moreover, the velocities are located 
where they are needed for the calculation of the convective fluxes of scalar flow variables. 

3.5 The Numerical Solution Procedure 
3.5.1 Numerical solution 

A control volume based fimte difference procedure [26] has been adopted in the 
present study to solve the governing equations numerically. The solution is initiated by 
descritizing the domain in a large number of non-overlapping control volumes. The governing 
equation is integrated around all such control volumes to yield a system of algebraic 
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Fig.Q.l) ; A typical three dimensional, scalar control volume and various 

nomenclature used in the cartesian co-ordinate system. 
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Fig.{3>3) : Scalar and staggered control volumes in three dimensional cartesian 
co-ordinate system. 
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equations, called the discretization equation [26]. The governing equations {Eqs.(3-l) to 
(3.8)) have both first and second order derivatives. As already mentioned, the first order 
derivatives are referred to the bulk convection term and the latter as the bulk diffusion term. 
These terms are then numerically integrated using the concept as applied to a combined 
convection-diffusion problem incorporating the hybrid differencing scheme as proposed by 
Patankar [26]. The pressure- velocity interlinkage in Eqs.Q.l) to (3.4) have been solved by an 
implicit finite difference procedure referred to as SIMPLE [26] (Semi-Implicit Method for 
pressure linked equations). A typical discretization equation for the general differential 
equation (£^.(3.14)) can be written as : 

in which, Ap , is the centre point coefficient ( see Fig.(3.\y) while Ap,Afp,Aj^,As,Ap and 
Ap are the neighbour point ( three dimensional ) coefficients, representing the combined 
effect of fluid convection and diffusion. Further, , represents the constant component of the 
linearized source term (viz., S = S^ + S^^^) while the slope has been embodied in the 
expression of Ap, the centre point coefficient. The centre point coefficient ^4^, is represented 
as ; 

A^=A^ *A^*A^*A^*Aj. -Ag-Sp.vol (3.16) 

Thus, in a system of n internal control volumes, n number of similar algebraic 
equations ( viz., Eq.Q.XS)) are obtained. Further, since the discretization equations are 
obtained from the same governing equation (£^.(3.14)), the former therefore embodies the 
same conservation principle as the latter one. This is an important feature of the control 
volume based formulation as against the routine Taylor series numerical solution procedure. 

Prior to the solution of the discretization equations the boundary conditions are also 
transformed into equivalent numerical form. Boundary conditions are only applied to control 
volumes located at the domain boundaries. Discretization equations are derived by the same 
procedure as mentioned. In essence, the incorporation of a boundary condition involves the 
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modification of and / or the source terms of the discretization equations of the boundary 
control volumes (see later). 

The resultant set of discretization equations for both boundary and internal nodes were 
solved using the Tri-Diagonal Matrix Algorithm (TDMA) adopting a line by line solution 
procedure. In this, a particular grid line, say in z-direction, is chosen and assuming the 
dependent variable (^ ) to be guessed in the x andy directions, the problem is reduced to a 
pseudo-one-dimensional situation and subsequently solved by TDMA. This was applied to all 
the grid lines in one direction and the entire process was repeated for the other two space 
directions to obtain a tentative solution of the dependent variable , (p . This typically 
constituted one iteration. The total number of iterations required was decided by the 
convergence criteria adopted, which in the present study was defined according to : 

E E E (317) 


in which. 


(3.18) 

The triple sum in Eq. (3. 17) represents the summation over the entire volume ( i.e., the 
calculation domain ). 

Under relaxation of the dependent variable was employed to achieve / enhance 
convergence. If is the value from the previous iteration of the dependent variable (pp {(pp 
implies (p evaluated at a point P )and a is the under relaxation factor then the change in the 
value of (pp can be slowed down from iteration to iteration by the following mathematical 
treatment [26] : 

<l>,-4>;^«( ^ -4>;) (3.19) 

Under relaxation factors for the problem chosen were determined by trial and error. Typically, 

CtMTRAL L'?RA 

A. 
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under relaxation factor of 0.5 was applied to u, v and w and while a value of 0.7 was used for 
the turbulence parameters, k and e. 

3.5.2 The computer program 

The concept of casting all governing differential equations to a general differential 
equation has been outlined in Section-3.3 . Consequently, a general purpose computer 
program was developed to solve £’^.(3.14). The computational procedure embodying the 
numerical solution outlined in the preceding section for the three dimensional model 
developed is shown in Fig.(3.4). The sequential steps of the computational procedure 
developed are as follows: 

i. The detailed information regarding the main grid, dimensions of the vessel or the 
section of the vessel under consideration, the fluid properties, turbulence model 
constants and the control parameters like convergence criteria and the maximum 
number of iterations are put in. 

ii. The grid arrangement for the scalar grid and the staggered grid is set up and the 
initialization of all variables ( u, v, w, k, p, p and e etc.) is done. 

iii. The w-momentum, v-momentum, w-momentum equations {Eqs. (3 .2) through (3 .4))are 
solved sequentially by solving the descritized equations for each of the control 
volumes by line by line solution technique. This is followed by solving the continuity 
equation (£^.(3.1)) embodying the SIMPLE algorithm [26]. The values of w, v and w 
for the calculation domain are corrected for the respective iteration. This is followed 
by solving the turbulence model equations (£^^.(3.5) and (3.6)) and the values of A: 
and € updated throughout the domain. 

iv. The fluid properties i.e. viscosity O^) are calculated (£^.(3.8))and updated for the 
calculation domain. 

V. Steps 3 and 4 are repeated till convergence criteria for all the flow, turbulence and the 
continuity equation is achieved. 

vi. Before the termination of the program the field distribution of flow and turbulence 
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Fig.(3A) : Flow chart for the three dimensional, turbulent flow calculation procedure. 
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parameters are printed. 


3.6 Results and Discussion 

The adequacy and appropriatness of the steady state, three dimensional, two phase 
turbulent flow model was assessed by carrying out some standard three dimensional flow 
simulation tests. This included flows in enclosed cavity with moving walls (F/g 5 ( 3 . 5 .a) and 
(3.5.b)), entrance length v. Reynolds number calculation in ducts of square cross-section and 
so on. These calculations were first carried out under laminar flow condition ( by setting // = 
//g and modifying and S„ appropriately in the flow equations). Subsequently, some 

of these tests were repeated under turbulent flow conditions ( using the k-e model). These are 
summarised below in detail. 

3.6.1 Laminar flow in an enclosed cubic cavity 

Steady, laminar , single phase flow in a three dimensional cavity can be conveniently 
represented by the Navier-Stokes equations presented in Section- 3.2. Thus, considering 
= n, or, =1.0 and S^= 5,,, = 0 in the calculation procedure, the computer program was 
run for a cubic cavity (/ = 0.5) filled with water {p, = 1000 kg/m^ and p = 0.001 kg/m.s). Two 
specific problems were solved for this configuration and this included : 

Prohlem-1: Flow with one of the cavity walls moving 

In this the south wall of the cubic cavity was made to move with a velocity of \m/s 
in the positive x direction (see Fig’.(3.5.a)). The remaining five solid walls namely the north, 
east, west, top and bottom walls were considered to be stationary. In order to define the 
boundary conditions for the given problem exactly, a zero velocity (no slip) condition was 
applied on the stationary walls. In addition to this, to incorporate the influence of the moving 
south wall on the fluid flow, the general descritization equation (£^.(3.15)) for the dependent 
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variable w, ( the x component of motion), i.e., 


(3.20) 

for the boundary control volumes lying in the immediate vicinity of the south wall were 
modified to : 

(3.21) 

where, S ^ = and u, is the given velocity of the moving wall( = 1 m/s). 

The predicted flow on the central vertical XY plane (at 2=0.25 m)is shown in 
Fig.(3.6). There as seen, the liquid adjacent to the moving wall accelerates along the direction 
of motion and turns upwards after hitting the east wall. Continuity requires that the flow 
along the north wall is in the -ve x direction and the fluid moves vertically down along the 
west wall. These form a recirculating loop with the eye of the vortex located at the geometric 
centre of the cavity. The predicted flow pattern as shown in Fig.{}.6) is clearly consistent with 
that one would expect for such a flow configuration. 

Problem-2: Flow with two parallel cavity walls moving in the same direction 

In this, both north and south walls of the cavity were made to move with a velocity 
of 0.01 m/s in the positive x direction (see Fzg.(3.5.b)). A similar treatment of boundary 
condition as outlined above was applied to u control volumes lying in the immediate 
neighbourhood of both the moving walls. It is important to mention here, that a substantially 
smaller velocity of moving walls was applied to this problem in contrast to the previous 
problem, primarily to assess the performance of the model developed for different order of 
magnitude of input variables. 

Figure (3.7.a) represents the fluid flow pattern on the central vertical plane at z = 0.25 
m . There as seen, the liquid is carried by the moving north and the south walls along with it. 
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Fig. (3. 7. Si) : The predicted flow pattern in the central XY plane (at z=0.25 m) with two 
parallel cavity walls moving.( Problem-!) 
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East Wall 





After hitting the wall ahead, the fluid flowing along the south wall turns upward whereas, the 
fluid flowing along the north wall turns downward. These two fluid streams as expected, 
meet exactly midway on the east wall (e.g., since both the north and the south walls of the 
cubic cavity are moving with the same velocity of 0.01 m/s in the same direction). After 
hitting each other the rising and falling fluids converge and turn backward at the centreline 
of the said plane towards the west wall. These as a consequence form two distinct contra- 
rotating, vortices one above and the other below the centre line . As seen the two loops appear 
to be exactly the mirror images of each other. Furthermore, the affect of the front and back 
walls on the fluid flow is practically non-existent as the central vertical XY plane is somewhat 
apart and equidistant from the front and the back walls. This essentially produces the two 
dimensional, symmetrical flow pattern. The affect of these solid walls on the fluid flow 
becomes noticeable when flow fields in two vertical, but identical planes close to front and 
back wall are examined. As seen in FigsQ.l.d) and (3.7.b) the two recirculating loops no 
longer remain symmetrical. Thus, in the vicinity of the front and back walls one of the loops 
gets restricted to a small pocket (comer region) in the main direction of flow. Once again the 
results for these two identical planes are seen to be the mirror images of each other. Predicted 
flow patterns shown in F/g.(3.7.a) through Fig.Q.l.c) are clearly consistent with that one 
would normally anticipate for such flow systems. 

3.6.2 Entrance length prediction for laminar flow 

When a flow enters a duct ( internal flow bounded by walls ) the viscous effects grow 
and boundary layer on the duct wall develops. The boundary layers eventually meet along the 
centreline and from this point onwards, the thickness of the boundary layers remain invariant 
in the streamwise direction. This is shown in Fig.{3.S.a) [33]. As seen, there is an entrance 
region where a nearly inviscid upstream flow converges and enters the tube. Viscous boundary 
layers grow downstream, retarding the axial flow u(x,r) at the wall and thereby, accelerating 
the centre core flow to maintain the incompressible, continuity requirement. At a finite 
distance from the entrance, the boundary layers merge and the inviscid core disappears. The 
tube flow then is entirely viscous, and the axial velocity adjusts slightly further until at x — L^, 
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Y 



North Wall 



Fig.QJ.h) : The predicted flow pattern in the XK plane (at z=0.05 m) with two parallel 
cavity walls moving. ( Problem-2) 
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North Wall 



Fig.(3 .7.C) ; The predicted flow pattern \n±Q XY plane (at 2=0.45 m) with two parallel 
cavity walls moving. ( Problem-2) 
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it no longer changes with x (no acceleration) and is said to be fully developed, i.e. u = u{r) 
only. Downstream of x = the velocity profile as well as the wall shear stress is constant, 
and the pressure drop varies linearly withx, for both laminar and turbulent flows (F 2 g.( 3 . 8 .a)). 

Dimensional analysis [33] shows that the Reynolds number is the only parameter 
affecting entrance length. Thus, 

L rt V n 

(3.22) 

For laminar flow [34], the accepted correlation for entrance length in ducts is ; 

L 

-^-KRe (3.23) 

Where A: is a constant of proportionality, the value of which, depends on the cross-section of 
the duct, e.g. for ducts with circular cross section, k = 0.06. In the present work entrance 
lengths were calculated for ducts with square cross-section (F/g.(3.8.b)) for laminar flow 
condition. 

The dimension of square cross-section (D) was varied between IxlO"^^ to 3x10"°^ m 
(dimensions in the y and the z direction ). Dimension of the duct in the x direction (the 
direction of main flow) was defined such that x »4. Zero velocity (no slip) boundary 
conditions were specified for north, south, top and the bottom walls (as already outlined in 
Section-3.2.3). Inflow and outflow boundary treatments were specified at the inlet (west) and 
outlet (east) planes respectively. Figure. {3. 9) represents a two dimensional grid arrangement 
near the entry and the exit planes. Since, u velocity (x component of motion) is known at the 
inlet plane, consequently, a numerical treatment analogous to that presented in Section-3.6 . 1 
can be applied to accommodate the exact inflow boundary condition into the calculation 
procedure. For all the grid points P lying adjacent to the outflow boundary, the coefHcient Ag 
is set to zero i.e. Peclet number is made sufficiently large { Pe " °°'> -^e "O) make the 
behaviour locally one way in this region. In addition to this, the average exit velocity is 
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Fjg.(3.8.a) : Schematic of the developing velocity profiles and pressure changes in the 
entrance of a duct flow. 
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Fzg.(3 .8.b) : Schematic of a duct with a square cross-section and the relevant co 
ordinate axis to the problem. 
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Inflow boundary 


• Scalar nodes 



Inflow boundary treatment 



Outflow boundary treatment 

Fig.(3.9) : The procedure adopted for modelling the inflow and outflow boundary 
condition for the prediction of entrance lengths. 
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specified at the outflow boundary on the basis of overall continuity which is updated from 
iteration to iteration by setting u{ni,j,k) = u(m-lj,k). 

Reynolds numbers (e.g., « inlet velocity ) were varied arbitrarily between 25 to 180 
and the flow field within the duct was predicted as a function of inlet Reynolds number. From 
the predicted results, entrance lengths were estimated graphically by plotting the centreline 
velocity as a function ofx. Entrance lengths thus estimated were non-dimensionalised (by 
dividing L, by D) and plotted against the inlet Reynolds number. This is shown in 10). 
As shown a straight line was obtained between the two, satisfying the requirement dictated 
by in Eq. (3. 23) (e.g., I D Re). This lends further confirmation to the appropriatness of 
the calculation procedure developed. 

3.6.3 Turbulent flow in an enclosed cubic cavity 

Follov/ing the tests under laminar flow situation summarised above, the three 
dimensional flow model was re-run under the turbulent flow condition for Problem-l ( flow 
with one of the cavity walls moving ). This was carried out by considering the exact form of 
£^5.(3. 1) through (3.4) in conjunction with the A:- e turbulence model ( Eqs.{3.5) through (3.8)) 
. Standard near wall treatment ( see Section-3.2.3) were provided for both flow variables and 
turbulence parameters. Numerical treatments on boundary conditions, at both moving and 
stationary walls were identical to those used earlier for Problem-l (see Section-3.6.1) 

It has been observed that the predicted flow in the central XY plane as well as 
elsewhere under turbulent flow conditions were qualitatively similar to those obtained under 
laminar flow conditions. The direction of flow recirculation, position of the vortex were 
practically identical for both laminar and turbulent flow conditions. To demonstrate the 
adequacy of the turbulent flow calculation results in Fig.(3.1 1), horizontal velocity u was 
plotted along the vertical central axis of the cubic cavity for a region close to the moving 
wall, for both laminar flow and turbulent flow predictions. There, as seen, for the region 
very close to the solid wall, the predicted velocity profiles are nearly identical under both 
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laminar and turbulent flow condition. This is because the shear stresses (at the wail) are 
expected to be identical for both the cases. This follows, since the shear stress is purely 
viscous and furthermore, that the viscous stresses are identical for both the situation, 
consequently, the velocity profiles have identical slopes there. However, in the region next 
to it, the momentum transfer normal to the moving wall is higher for the turbulent flow 
because of the presence of additional stresses i.e. Reynolds stresses. This makes the turbulent 
flow profile less steeper than the corresponding laminar flow situation. The velocity profiles 
in a nearly similar situation i.e. stationary wall with the bulk liquid moving with a velocity of 
1 m/s have been studied [35] and the corresponding results for the near wall region is 
reproduced in F/g.(3.12). Qualitative agreement between the present estimates and literature 
is at once evident. This lends further credence to the three dimensional turbulent flow 
calculation procedure developed in this study and consequently, the model can be now 
extrapolated with confidence to investigate melt flow behaviour in metallurgical tundish 
systems. 
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CHAPTER -4 


MODELLING OF FLOW PHENOMENA IN RECTANGULAR 

SHAPED TUNDISHES 


In the previous chapter, the salient features of upgradation of a two dimensional 
calculation procedure [ 1 8 ] to a steady state, three dimensional, two phase turbulent flow model 
and its validation with respect to some benchmark solutions were presented. In this chapter the 
three dimensional turbulent flow model is applied to investigate fluid flow and residence time 
distribution (RTD) in rectangular shaped tundishes with single entry and exit and no flow 

control devices. 


4.1 Introduction 

A continuous casting set up essentiaUy consists of three components, the teeming ladle, 
the m oHieh and the continuous casting mould. The tondish is the intetmediate vessel between 
the ladle and the mould which fimctions as a distributor and a buffer vessel and helps mamtam 
a relatively uniform flow as the teeming ladle is being emptied and replaced by a new one. 
f/gare.(4.1) shows a schematic of the continuous casting process. Over the years the quality 
requirements of the final steel product have become very stringent and the tundish offers the 
last opportunity for the removal and modification of impurities and the adjustment of tnm 
cWicI composition by micro-alloying. With the emphasis on the novel casting methods like 
thin slab and strip casting in a single step, the role of tundish has become even more cntic 

The performance of a tundish is very closely related to the nature of steel flow i.e. a well 
mixed flow would be more suitable for dissolution and alloying where as a plug flow is more 
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F/g.(4.1) : Schematic representation of the continuous casting process. 
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suited for inclusion separation. Therefore, tundish design is important for the tundish to 
perform the desired functions effectively. Towards this, mathematical modelling can give a 
quantitative assessment of flow and turbulence parameters and the associated RTD distribution 
and thus provide useful insight into the overall process. This intum can help achieving an 
optimal tundish design. 


4.2 Characteristics of Metallurgical Tundish System 

When fluid flows through the tundish it is important to establish the time spent in the 
system by individual fluid elements. The average time or the nominal holding time, of the fluid 
in the system is easily calculated from the definition 

Volume of fluid in vessel 

(4.1) 

Volumetric rate of fluid flow 

However, some fluid individual elements may spend a longer, and others a shorter 
period of time in the system. The minimum amoimt of time a fluid spends inside the vessel is 
called the breakthrough time. This departure of actual residence times from the mean i.e., the 
distribution of residence times is an important characteristic of the tundish system and 
influences the latter's performance as a reactor. The residence time distribution of a fluid 
flowing through the tundish can be determined by means of tracer studies. This involves the 
addition of the tracer to the stream entering the vessel, and then measurement of its 
concentration at the exit with respect to time. 

Several methods have been developed for introducing the tracer material into the vessel, 
but the two most important [25,36] are : 

i. continuous addition of tracer to the input stream, starting at a particular time (step 
input), and 

ii. addition of the tracer over a short time interval, the duration of which is negligible in 
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comparison with the mean residence time of fluid in the vessel (pulse input). 


Of the above two methods, pulse input of tracer is more popular as in many tests the use of step 
input of tracer may lead to extensive contamination of product or alter the fluid flow behaviour 
itself. It has also been shown [36] that the information derived from pulse addition of tracer is 
ideal for discriminating the hydrodynamics of different tundish systems. 


The results obtained from pulse addition of tracer are plotted in the form of 
dimensionless concentration vs. dimensionless time and are called C-curves. 




c 

Qiv 


(4.2) 


where, 

Q = quantity of the tracer injected, 

V = volume of the fluid in the vessel. 

The denominator of Eq.{A.2) represents the average concentration that the tracer would 
reach if it were perfectly mixed within the vessel. The area under the C-curve is umty as all 
tracer entering the system should eventually leave it ( f Cddd =1). From the C-curve 
parameters like the minimum breakthrough time ( ), the time at which the peak 

concentration is reached ( average residence time ( /^) can be obtained . These 
can be converted into their equivalent dimensionless forms ( = /„,„/ , 6^^ 

= t^g/ o where t, is the theoretical residence time as defined by Eq.{A.\). To interpret the 
results from the C-curve quantitatively the vessel volume can be assumed to be made up of 
three parts namely, a well mixed volume (V^, a dispersed plug volume (F^) and a dead volume 
(VJ [36]. The fiactional volumes can be calculated as [36] : 

F, = (l-9^ (4.3) 
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(4.4) 



The relative amounts of these volumes determine the efficiency of a given tundish and can also 
be used as appropriate criteria for improving the tundish design or for comparative analysis. 
The preceding formulation assumes that only negligible volume of liquid spends a time greater 
than twice the theoretical residence time in the tundish. 


4.3 Mathematical Modelling of Tundish Hydrodynamics 

As mentioned already, a metallurgical tundish may be regarded as a trough into which 
molten metal is being poured from a teeming ladle at one location and discharged from the 
tundish through one or more outlets. The flow is known to be highly turbulent in the vicinity 
of the pouring stream while elsewhere in the tundish, flow is in the transitional regime [36]. 
Further, the nature of the fluid flow is essentially three dimensional and elliptic. Except for the 
initial transient and the final drain out stages, the flow induced in a tundish ( when bath height, 
exit velocity etc. remain constant ) can be described mathematically by the steady state, three 
dimensional calculation procedure outlined in Chapter-3. Once flow variables and turbulence 
parameters are predicted, these can be coupled with an appropriate convection-diffusion 
equation to predict the RTD (residence time distribution) pattern in the tundish. (see later) 

In the present study, fluid flow behaviour in a model tundish used by Singh and Koria 
[17] with single inlet and outlet and no flow control devices has been analysed. The 
characteristic design features of the model tundish [17] have been summarised in Table-A.\. 
To mathematically model the tundish the sloping side walls were neglected and the tundish as 
a first approximation was assumed to be perfectly rectangular. The schematic of the tundish 
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modelled [17] in the present study has been shown in Fig.{A.2). 


Table-4. 1 : Characteristic parameters of the single strand model tundish [17]. 


Parameter 

Model Tundish 

Fluid 

Water, (//=0.001 kgm'^s'\p=l000 kg/m^ 

Base length, m 

1.0 

Width (at base), m 

0.30 

Height, m 

0.37 

Liquid depth, m 

0.26 

Shape 

rectangular with sloping walls (15°) 

No. of strands 

one 

Ladle shroud diameter, m 

0.021 

Volumetric flow rate, mVs 

1.55 X 10-^ 

Inlet Reynolds number 

9.0 X 10^ 


4.3.1 Modelling assumptions 

The following assumptions and idealisations were made in the present work to 
investigate the fluid flow behaviour in tundish. 

i. The flow is assumed to be turbulent and macroscopically steady i.e., the effect of 
transience during filling and emptying of tundish was neglected. 

ii. The melt surface is assumed to be flat and phenomena like wave formation, presence 
of slag layer at the free surface and vortexing during drainage etc. are neglected. 

iii. Entrainment of air or gas through the metal inlet is neglected (the ladle stream is 
essentially shrouded ). 
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iv. The flow is essentially Newtonian and incompressible and 

V. The system is isothermal and thus the effect of temperature differential on the melt flow 
is neglected. 

4.3.2 Governing equations and boundary conditions 

As already mentioned, the nature of fluid flow in tundishes is three dimensional and 
can be represented by the continuity equation and the three components of the Navier-Stokes 
equations already outlined in Chapter-3 (Refer £'95.(3.1) to (3.4)). The effective viscosity 
appearing in these equations can be obtained from the model of turbulence (Refer Eqs.(3.5) 
to (3.8)). The applicable set of boundary conditions to the tundish system are : 

i. At the free surface ( in the north direction ) zero shear is assumed to be transmitted for 
the velocity components. Zero gradient of k and e were specified at the free surface. In 
conjunction with these, vertical component of flow is set to zero. Mathematically, 

Aty = 0.26 m, 0sx<l.0m and 0 i z ^ 0.30 m, 

v = 0, di/<^ = <9w/c>y = 0 , = (9€/(^ = 0 

ii. No slip condition for the velocity components was specified at the boundary walls, (see 
Fig. (4.2)) The values of k and e were set at zero at the wall. 

Aty = 0 , Osxsl.Om and 0 s z ^ 0.30 m, 
u = v = w = 0 , fc=e = 0. 


At X = 0 and x = I m. Os ys 0.26 m and 0 s z s 0.30 m, 
u = v = w = 0 , k=e = 0. 


Atz = 0 andz = 0.3m. Os ys 0.26 m and Os xs 1.0 m, 
M = v = w = 0 , k—e. — 0. 
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In addition to this, at the exit, outflow boundary condition was specified by considering 
relevant Peclet number to be mfinity (e.g., zero diffusion). Regions close to the wall have been 
modelled using the standard wall function procedure [21,32] (for details see Chapter-3) In 
addition to these, at the inlet the normal velocity was specified from the total volume flow rate 
( v=0.41 m/s). The inlet values of k and 6 were estimated from the following relations [37]. 

= 0.01m; (4.6) 

. 1 - ' /A 


where u, is the normal inlet velocity and d, is the diameter of the inlet nozzle. At the outlet 
nozzle, the normal velocity was calculated from the overall mass balance [38]. 

Computations were performed usir^ uniform grid covering the whole tundish. The grid 
set of 52 X 12 X 7 were chosen for the x (longitudinzil), y (vertical) and the z (transverse) 
directions respectively. Under-relaxation factors of 0.5 for velocities and 0.7 for scalar 
quantities were adopted ( see Chapter - 3 for details on the under-relaxation procedure 
employed ). A typical execution time required about 150 min of CPU time on a /7P-9000 
machine. 


4.3.3 Governing equation of tracer dispersion 


As pointed out previously, the residence time distribution (RTD) characteristic of a 
tundish can be studied effectively by tracer dispersion studies. Thus, following a pulse addition 
the appropriate conservation equation of the added tracer can be described in terms of the 
following convection - diffusion equation in cartesian co-ordinate system: 


|-(Q . ^(uC) . j-(yC) . |-(wQ - * j-(D,^) ^ 

dt dx dy dz dx dx dy dy dz dz 


(4.8) 
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diffusivity and Dj is the eddy diffUsivity. The latter, can be directly equated to the turbulence 
kinematic viscosity assuming turbulent Schmidt number to be approximately unity 
i.e. fij! (pDj) = 1 (valid for a wide range of engineering problems). A zero flux boundary 
condition is imposed on all the solid surfaces and the free surface. At the exit, for the control 
volume passing through it, a sink term was incorporated into £’^'.(4.8) as: 

dC 

— ^ * div (uC) = div (D grad C) - KC_ (4.9) 

dt 

where K = { cross-sectional area of the outlet face of control volume) x (the normal velocity 
through the outlet face ). The negative source (e.g., sink) term, -KC, was dumped in (i.e., 
Ap = YAnb ~ numerical stability. In addition to this, an initial condition, i.e., 

atr = 0,Q = C; atI=S,J=lO,K = 4 

was applied to Eq.(4.S). The specified I, J and K coordinates are respectively the location of 
the tracer input, where, C° = 1000 kg /m^ (tracer having the same density as the bulk liquid) 
was used in the numerical solution scheme. It is instructive to note here that one way coupling 
between £^.(4.8) and the flow and turbulence model equations allows flow computation to 
precede the calculation of the tracer dispersion. Once converged flow field is obtained, this was 
applied to predict the C, ( x,y,z,t) field. A convergence criteria of lO"®^ was qiplied on £^^-.(4.8) 
and (4.9). In addition for every time step a minimum number of iterations ( = 7) was carried out. 
From the predicted concentration at the outlet as a function of time, t, the C-curve can be 
conveniently constructed as explained in Section-4.2. 


4.4 Results and Discussion 
4.4.1 Fluid flow behaviour 

The predicted velocity field in the rectangular shaped tundish has been shown as a 
series of 2D plots (e.g., vertical planes and horizontal XZ planes ) in Figs(4.3) and (4.4). 
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As seen from these figures, the general nature of the flow is of three dimensional nature with 
pronounced spatial variations. 

Figure (4.3 .a) represents the flow on the central vertical XY plane passing through the 
inlet and the outlet. There as seen, the incoming jet of water flows downward, entrains liquid 
from the surrounding and thus decelerates considerably as it proceeds towards the base of the 
vessel. The flow in and around the inlet region is highly turbulent (predicted is of the order 
of 0.12 kg/m.s). The fluid after hitting the base, moves outwards towards the vessel side wails. 
Closer to the left wall, this outwardly flowing fluid rises to the free surface and finally, turns 
inward towards the jet entry region. This intum forms a small but distinct recirculating loop 
(see the left portion of Fig.(4.3.a)). To the right of the inlet stream, the fluid in the upper region 
of the vessel is seen to be converging towards the jet entry, while in the vicinity of the tundish 
base, the flow is away from the jet impinging region. It is also clear from Fig.{4.3 .a) that the 
velocities drop significantly with increasing distance from the inlet region. Fluid velocities 
again become appreciable near the outlet. 

The fluid.after hitting the vertical walls (those along the z direction) also moves upwards 
towards the free surface as seen from Fzg.(4.3.b). There the velocity field has been plotted on 
a vertical A7 plane in the immediate vicinity of the wail. As seen, close to the inlet region, the 
flow is directed vertically upward ( a consequence of the continuity requirement) while in the 
vicinity of the exit, the flow is directed vertically downwards. As one might anticipate, the inlet 
and the outlet in the tundish strongly affect the flow pattern in their immediate vicinity. 


Similarly, when results are plotted on horizontal planes (see F/g.(4.4.)), it is readily 
seen that the inlet jet after hitting the base of the tundish spreads in all directions, but 
predominantly towards the front and rear walls ( walls along the z direction ). Furthermore, the 
jet spreads very little towards the exit nozzle (F/g.(4.4.a)), as it is much more easier for the 
liquid to move towards the side walls (e.g., the momentum of the liquid is much higher in this 
direction). In the central region of the tundish the flow is noticeably weak. Finally, in the 
proximity of the exit nozzle the flow becomes strong again because of the drainout action. In 
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0.07 m/s 


FigiA3 .a) : The predicted flow pattern in the central vertical XY plane (at z = 0. 1 5 m). 
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/’/^.(4.3.b) : The predicted flow pattern in a vertical JO^plane close to the wall f at z = 
0.03 m) . 
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/r/g.(4.4.a) : The predicted flow pattern in a horizontal XZ plane close to the base of the 
tundish (at>^ = 0.012 m). 
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addition to these, weak recirculation of the liquid is seen to the right of the exit point (e.g., 
between the exit nozzle and the wall to its right). 

s 

In Fig.(4.4.h), fluid flow pattern has been plotted on the horizontal XZ plane passing 
through the middle of the tundish ( at>^.l 1 m). There the liquid can be seen to be converging 
very strongly towards the entry nozzle particularly from the front and the rear of the inlet. The 
flow vectors in the central region of the plane are directed towards the inlet region whereas, 
fluid flow closer to the walls is directed towards the outlet . A weak recirculating loop is also 
seen to be forming about the outlet. Finally, in F/g.(4.4.c), fluid flow pattern is shown on the 
horizontal XZ plane (aty = 0.24 m) close to the free surface of the liquid. There also the liquid 
is seen converging towards the entry jet region from left part of the vessel. Immediately, above 
the outlet region, the flow is practically absent in this plane (free surface is farthest from the 
exit). 


From the preceding discussion, the following general observations can be made about 
the nature of the fluid flow in rectangular shaped tundishes with single entry and exit nozzles 
and no flow control devices. The entering jet spreads outwardly following impingement and. 
this as a consequence, entrains the surrounding liquid. Subsequently, the liquid flows towards 
the walls of the tundish and then, towards the free surface and finally converges to the 
entering jet of water. This characteristic of flow is identical for all directions in the vicinity of 
the inlet nozzle only. The flow pattern however, gets considerably modified in the vicinity of 
the exit, where the flow direction is seen to be governed predominantly by the drainout action. 
These features of flow in metallurgical tundish systems which is the result of the combined 
action of an impinging jet and a outgoing stream are now rather well established. Present flow 
predictions are qualitatively in line with similar computational studies reported earlier for the 
tundish system. To substantiate this, computational results of Sahai from ref.39 is shown in 
Fig.(4.5). 
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Fig.{AA.h) ; The predicted flow pattern in a horizontal ^ plane (at>^ = 0.11 m). 
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0.02 m/s 


Fig.{AA.c) : The predicted flow pattern in a horizontal XZ plane close to the free 
surface. {zXy - 0.24 ni). 
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(a) near the plane of symmetry and (d) closer 
to the inclined wall . 


Fig.i4.5) : Predicted flow field in longitudinal, vertical planes with no flow control 

devices in a tundish with inclined wall [391. 


4.4.2 Residence time distribution characteristics 


As pointed out already, the residence time distribution of the fluid in the model tundish 
was investigated by pulse addition of a tracer at the fluid inlet location i.e. at t = 0 , the inlet 
control volume was made full of tracer of the same density as that of bulk liquid phase (e.g., 
water). The concentration of the tracer was then computationally monitored as a function of 
time at the control volume passing through the exit nozzle and one control volume vertically 
above it. The concentration of tracer at these locations was subsequently non-dimensionalised 
by dividing it by the average concentration that the tracer would assume in the vessel, if it were 
perfectly mixed in the absence of any drainout ( Q = {^mcJ^water) ~ kg/m^). Similarly, the 

time was non-dimensionalised using the nominal holding time (see £'^.(4.1)) of the fluid in 
the vessel ( t^). The result thus obtained ( C/Co vs. t! was then plotted to get the C-curves. 
In Fig.iA.6), the C-curves obtained computationally for the two control volumes (1 mm and 3 
mm above the exit nozzle) are shown. There as seen, the tracer is initially of negligible 
concentration at the exit stream, but it quickly increases to a maximum, before showing a 
gradual fade out. This characteristic of the RTD curve is similar to those obtained earlier by 
various investigators for the tundish systems [17,36,39]. The breakthrough time was calculated 
as the time when the concentration in the vicinity of the exit reached 0.1 % of the average 
concentration of the tracer (i.e., C = 4 x kg/m ^ ). 

Numerically predicted breakthrough time and the time at which peak concentration 
is achieved at the outlet together with those obtained by Singh and Koria [17] have been 
summarised in Table-4.2. As seen there, considerable differences appear between observed and 
predicted peak concentration time. 

The possible reasons for the discrepancy between the observed and predicted peak 
concentration time can be attributed to: 


sensitivity of the results to the step height df and the numerical grid chosen, 
inadequacy of the model of turbulence to accurately simulate such systems and 
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Fig.iA,6) ; Predicted RTD curves for monitoring location at \mm and 3 mm above the 
nozzle exit. 




finally, 

iii. the assumption of a perfectly rectangular shaped tundish (e.g., ignoring the wall taper). 

Table-4.2, \ Comparative results of RTD characteristics obtained via physical modelling [17] 
and mathematical modelling. 


RTD Characteristics 

Mathematical Modelling 

Physical Modelling [17] 

rngm 

3mm® 

Breakthrough time 

25 

20 

30 s 

Peak concentration time 

445 

465 

260 s 


® vertical distance of monitoring location above the exit nozzle. 


The time dependent tracer dispersion equation ( £^.(4.8)) has been solved by the fully 
implicit marching integration procedure [26], which is known to be relatively more accurate for 
larger time steps. Hence, the numerical calculations were carried out considering two time 
steps (5 and 10 5 ) to ascertain the influence of step height on the predicted results. As seen in 
Fig.(4.7.a), increasing the time step two fold, from 5 to 10 s, does not bring in any change in 
the solution. Consequently, results shown in Table-4.2 can be assumed to be independent of At 
chosen. Similarly, two different grid systems 25 x 8 x 7 and 52 x 12 x 7 were used. Only 
marginal changes in predicted flows and concentration fields were observed. 

The fluid flow in a steelmaking tundish is known to be turbulent only in the inlet 
region [36]. Further, the intensity of flow variables and turbulence parameters in tundishes are 
essentially determined by the kinetic energy of the incoming jet and the vessel geometry. The 
flow is in the transitional regime at locations away from the inlet stream. To assess the role of 
a turbulence model in the calculation procedure, the bulk effective viscosity model proposed 
by Pun and Spalding [40] was applied in the numerical calculation procedure in leu of the k-e 
model of turbulence. Pun and Spalding proposed the following empirical correlation, through 
dimensional reasoning for computing turbulent flows in sudden expansion geometry. 
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(4.10) 


In £^.(4.8), jUg , the dependent variable is related to the medium density p, the reactor 
length H, the reactor diameter D and the kinetic energy rate of incoming fluid, {mu^. is an 
empirically fitted constant and assumes a value 0.012 for sudden expansion geometry. The 
calculated via the above correlation ( assuming, D = width along the z direction and H= depth 
of the liquid) for the present study was found to be approximately 0.04 kg . In contrast, 

the k-e model of turbulence predicted in the range of 0.05 to 0.2 kg Thus, 

calculations were repeated using the bulk effective viscosity model [40] whereby, p^ was 
arbitrarily specified in the range of 0.05 to 0.1 kg ni^s'‘ within the tundish. The affect of 
effective viscosity on the residence time distribution is shown in Fig.{A.lh). There as seen 
numerical value of p^ shifts the RTD curve and consequently, the residence time of fluid in the 
tundish appreciably. However, RTD characteristics like breakthrough time and peak 
concentration time remain virtually unaffected. Hence, the choice of the turbulence model does 
not appreciably affect the predictions regarding breakthrough time and the peak concentration 
time. Evidently, the total residence time which is a measure of the area under the curve does 
change significantly with different choice of turbulence model. 

Finally, as mentioned already the wall taper of 15® has been neglected in the present 
niunerical simulation. This has lead to an underestimation of the tundish volume by over 1 8 %. 
Hence, the breakthrough time predicted in the present study is lower fi'om that predicted by the 
physical modelling route [17]. Similarly, the poor agreement between the predicted and 
observed peak concentration time appears to be related to the fact that tapering of sidewalls 
alters the flow significantly, and makes it predominantly directed along the base of the timdish 
[17,39] which has not been computationally observed in the present study. Quantitative 
comparison thus is meaningful only when the exact vessel geometry is incorporated into the 
numerical solution scheme. This is however left for a future work. 
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Fig.{A.l.h) ; The affect of numerical value of effective viscosity on the prediction of 
RTD characteristics in the tundish system. 


94 


CHAPTER -5 


CONCLUDING REMARKS 


The following three major conclusions can be derived from the present study: 

i. A hydrodynamic model has been developed, by using a two phase, two dimensional 
turbulent flow model in conjunction with an appropriate heat transfer correlation to 
estimate the melting rates of solids in gas stirred metallurgical reactors. It was 
demonstrated that the mathematical model simulates experimental observations 

^ realistically. The predictive capabilities of the model were also evaluated with 
references to various process and numerical parameters. It has been shown explicitly 
that for accurate prediction of rate processes in such two phase systems, realistic inputs 
of bubble diameter, drag coefficient values etc. are required. 

ii. In the absence of computational facility, a simple, yet effective macroscopic equation 
has been suggested for estimating melting rates of solids in the plume region of gas 
stirred baths. It has been demonstrated that the macroscopic model produces estimates 
that are comparable to the detailed differential model results. 

iii. Finally, the two dimensional, two phase calculation procedure was modified and a 
steady state, two phase, three dimensional, turbulent flow model has been developed, 
based on the control volume formulation, embodying the SIMPLE algorithm. The three 
dimensional calculation procedure developed has been tested against standard test 
solutions and qualitative agreement shown. The model was shown to simulate 
qualitatively melt flow behaviour and residence time distribution in rectangular shaped 
tundish systems. 
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CHAPTER -6 


RECOMMENDATIONS FOR FUTURE WORK 


The scope of the present work can be further expanded considering the following: 

i. In the present study, fluid flow has been evaluated in rectangular shaped tundishes 
with no flow control devices . The three dimensional turbulent flow model can be 
conveniently modified to study transport phenomena in tundishes with inclined walls 
and flow control devices (e.g., dams, weirs, slotted baffles etc.) 

ii. Argon and/or nitrogen injection which is increasingly being applied in tundish 
operations is a relatively new area that can be investigated with the aid of the present 
model. No comprehensive two phase simulation (either experimental or theoretical) of 
metallurgical tundish system has been reported in the literature to-date. 

iii. The two phase, three dimensional, turbulent flow model which has been developed in 
the present work can be suitably modified, to carry out equivalent computations in 
polar co-ordinate system. The code can then be used to simulate a variety of 
metallurgical processes such as, asymmetric gas injection in vessels, side blown 
converters, -degassers and so on. 

iv. ' The influence of temperature variations on fluid flow (firee convection affects) has been 

neglected in the present study. Integration of heat transfer aspects to the present model 
will enable it to predict the affect of temperature losses or auxiliary heating on fluid 
flow. Similarly, appropriate form of Maxwell equations can be incorporated in the body 
force terms of the Navier-Stokes equation in the present model, to study the affect of 
magnetic field on fluid flow in metallurgical tundishes. 
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